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Published  15  February  2016  Characterizing  the  fast  evolution  of  microstructural  defects  is  key  to  understanding  "crackling" 
phenomena  during  the  deformation  of  solid  materials.  Forexample#  it  has  been  proposed  using 
atomistic  simulations  of  crack  propagation  in  elastic  materials  that  the  formation  of  a  nonlinear 
hyperelastic  or  plastic  zone  around  moving  crack  tips  controls  crack  velocity.  To  date,  progress  in 
understanding  the  physics  of  this  critical  zone  has  been  limited  due  to  the  lack  of  data  describing  the 
complex  physical  processes  that  operate  near  microscopic  crack  tips.  We  show,  by  analyzing  many 
acoustic  emission  events  during  rock  deformation  experiments,  that  the  signature  of  this  nonlinear 
zone  maps  directly  to  crackling  noises.  In  particular,  we  characterize  a  weakening  zone  that  forms 
near  the  moving  crack  tips  using  functional  networks,  and  we  determine  the  scaling  law  between  the 
formation  of  damages  (defects)  and  the  traversal  rate  across  the  critical  point  of  transition.  Moreover, 
we  show  that  the  correlation  length  near  the  transition  remains  effectively  frozen.  This  is  the  main 
underlying  hypothesis  behind  the  Kibble-Zurek  mechanism  (KZM)  and  the  obtained  power-law  scaling 
verifies  the  main  prediction  of  KZM. 

The  spatio-temporal  evolution  of  acoustic  signals,  a  class  of  crackling  noise1,  is  a  direct  result  of  failing  atomic 
bonds  during  material  fracture.  Such  signals,  if  properly  interpreted,  may  be  used  to  better  understand  the 
dynamics  of  rupture  progress  in  the  vicinity  of  crack  tips  over  a  broad  range  of  scales  and  conditions2-4.  During 
microcrack  propagation,  part  of  the  stored  energy  near  the  crack  tip  is  consumed  in  the  breaking  of  molecular 
and  atomic  bonds,  resulting  in  new  crack  surface  area.  The  key  to  understanding  the  crackling  process  lies  in  the 
characterizing  structure  of  the  near- tip  region  of  such  micro  cracks,  where  stress  amplitudes  are  large.  Due  to  the 
microscopic  size  and  high  speeds  encountered  in  the  vicinity  of  the  crack  tip,  direct  measurements  are  difficult, 
and  analysis  typically  relies  on  computational  techniques5.  Because  of  the  large  strains  present  near  crack  tips, 
nonlinear  elastic  and/or  inelastic  contributions  must  occur,  and  recent  work2,5  suggests  that  non-linearity  around 
the  moving  crack  tip  governs  the  rupture  velocity  Specifically,  the  local  hyperelastic  or  plastic  zone  around  the 
moving  crack  tip  enhances  energy  flow  for  stiffening  systems,  and  reduces  energy  flow  for  softening  systems, 
resulting  in  increases  and  decreases  in  the  fracture  velocity  respectively2.  Furthermore,  investigations  of  “s/ow” 
cracks  in  gels  demonstrate  a  link  between  the  spatial  energy  flow  around  the  rupture  tip  and  the  curvature  of  the 
tip5.  This  link  is  thought  to  be  responsible  for  inaccuracies  in  linear  elastic  analyses  that  are  commonly  used  in 
material  science  for  simulating  crack  tip  processes,  with  applications  ranging  from  metal  fatigue  to  earthquake 
nucleation.  In  all  of  the  mentioned  studies  (but  see6),  the  assumption  is  that  cracks  evolve  under  equilibrium 
conditions,  and  all  crack  tip  processes  satisfy  an  (quasi)adiabatic-equilibrium  assumption.  Under  this  assump¬ 
tion,  the  ramp  time  (the  duration  that  the  system  approaches  the  critical  point  at  which  unstable  crack  growth 
proceeds)  is  at  least  of  the  order  of  the  relaxation  time  of  the  system7-10.  In  this  study  we  question  this  assumption 
by  showing  that  fast  processes  (including  many  sub -excitations  possibly  due  to  lattice  distortions  or  dislocations) 
involved  in  the  structure  of  fast-moving  crack  tips  induce  a  “frozen”  regime  where  -  for  a  very  short  time  -  the 
crack  tip  is  out  of  equilibrium  and  the  dynamics  are  very  slow. 

In  this  study  we  use  a  functional  network  method11,12  applied  to  acoustic  emission  (AE)  data  recorded  during 
rock  deformation  and  rock  friction  experiments  (see  methods  section)  to  show  that  moving  microcracks  contain 
signatures  of  non-linearity  We  discover  that  the  onset  of  the  non-linear  stage  prior  to  unstable  failure  coincides 
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Figure  1.  Q-profiles  representing  dynamic  crackling  noises,  (a)  Three  main  stages  of  typical  acoustic 
crackling  noises  as  shown  in  normalized  Q-profiles:  S-W-D  phases  correspond  with  strengthening,  weakening 
and  decelerating  stages,  respectively,  (b)  A  typical  -R  profile  as  generated  from  dataset  Lab.EQl.  (c)  A  schematic 
representation  of  impulse  zone  in  transition  from  S  to  W  phase,  (d)  Five  recorded  stick-slip  events  with 
(dynamic)  strain  gauges  measurements  in  centimeter  scales  rock- interfaces12,29,  (e)  An  event  from  (d)  where 
stresses  dynamically  drop  about  25  MPa  during  a  stick-slip  experiment. 


with  the  nucleation  of  kink  instabilities,  and  we  use  “network  strings”  to  visualize  spatio-temporal  evolution  of 
these  topological  defects.  For  the  first  time,  we  show  that  emitted  crackling  noises  hold  the  signature  of  the 
Kibble-Zurek  mechanism  (KZM)7"10  that  provides  an  estimation  of  the  defect  density  as  a  function  of  the  traversal 
rate  across  a  phase  transition.  In  the  current  case,  the  “phase  transition”  is  a  transition  between  a  strengthening 
and  unstable  weakening  state  during  microcrack  propagation.  A  key  output  of  KZM  relates  faster  ramp  rates  to 
higher  defect  density,  as  the  result  of  a  spontaneous  symmetry  breaking  process.  We  measure  real-time  evolution 
of  a  first-order  correlation  function  of  the  system  (in  network  space)  and  verify  the  main  prediction  of  KZM: 
namely  the  power-law  scaling  of  the  (frozen)  correlation  length  with  the  ramp  rate.  Moreover,  we  show  that  the 
correlation  length  near  the  transition  remains  effectively  frozen.  This  is  the  main  underlying  hypothesis  behind 
the  Kibble-Zurek  mechanism7,8.  In  addition,  using  our  laboratory  datasets  we  illustrate  that  the  adiabatic-impulse 
transition,  as  the  core  of  the  KZM  hypothesis,  can  be  used  to  infer  an  approximate  weakening  rate. 

We  analyse  AE  waveforms  (the  laboratory  analogue  of  seismograms  due  to  rock  fracture  or  earthquake 
rupture)  under  different  simulated  depth  (pressure)  conditions  and  loading  paths,  and  on  two  rock  types: 
Westerly  granite  and  Basalt  from  Mount  Etna,  Italy  (datasets  Lab.EQl,  Lab.EQ2,  Lab.EQ3  and  Lab.EQ4  -  see 
methods  section  and  supporting  information).  We  apply  tools  from  the  theory  of  complex  networks  to  ana¬ 
lyze  emitted  noises  from  microscopic  cracks,  where  the  acoustic  time  series  recorded  at  each  sensor  is  repre¬ 
sented  as  a  node  [Si-section  1;  refs  1 1,12].  These  results  allow  us  to  develop  an  interpretation  of  recorded  multiple 
acoustic- crackling  signals  involving  a  microsecond  evolution  of  different  dynamic  crack  tip  phases  as  encoded  in 
the  network  modularity  (Q)  which  we  refer  to  as  “Q-profiles”  (see  methods  section).  This  evolution  can  be  broken 
down  into  three  distinct  phases  (Fig.  la,d,e):  (1).  The  S-phase:  an  initial  strengthening  phase  preceding  the  critical 
point  at  which  point  weakening  and  catastrophic  failure  begins;  (2)  the  W-phase:  a  fast-slip  or  weakening  phase; 
and  (3)  the  D-phase:  a  slow  slip  or  decelerating  phase11-13.  To  better  understand  the  S-phase  and  how  the  transi¬ 
tion  occurs  across  the  critical  point  between  S  to  W,  we  use  the  reciprocal  of  modularity  (R=  Q_1)  profiles  (i.e., 
“R-profiles”)  which  closely  resemble  dynamic  stress  profiles  commonly  used  to  characterize  rock  failure  (Fig.  1). 
R-profiles  are  indicative  of  the  dynamic  stress  changes  due  to  a  given  cracking  event.  Rmax  corresponds  to  the 
critical  point,  where  the  failure  occurs  and  fast-weakening  begins. 

In  the  following  discussion,  we  describe  the  observation  of  “defect”  formation  prior  to  onset  of  the  W  phase. 
To  proceed,  we  define  the  critical  zone  onset,  Rc,  as  the  value  of  R  at  the  time  of  the  first  impulse  in  the  inverse  of 
mean  betweenness  centrality  (1  /log  ( B .  C.) )  profile,  where  ()  indicates  the  mean  value  of  all  nodes  (Figs  2  and  3). 
In  Fig.  2,  we  show  1  /log  ( B .  C. )  profiles  for  6  acoustic  events  from  our  laboratory  tests.  For  all  events  this  profile 
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Figure  2.  (a,b)  Two  typical  acoustic  emission  events  from  cracking  Granite  samples  (from  Lab.EQ3).  We  have 
shown  scaled  recorded  acoustic  waveforms  in  -800  ps  and  the  corresponding  normalized  Q(t).  (c)  Different 
events  from  our  tests  with  the  signature  of  inverse  of  mean  betweeness  centrality  which  shows  divergence  of 
the  parameter  in  vicinity  of  the  nucleation  zone.  Based  on  the  resolution  of  our  measurements,  the  total  time  of 
order- disorder- order  transition  sequences  stretches  between  0.5-16  ps 


is  characterized  by  a  narrow  interval  at  the  transition  from  the  S  to  W  phase.  Later,  we  will  show  that  the  first 
impulse  corresponds  to  the  first  nucleated  defect  in  transition  from  S  to  W  while  the  second  spike  indicates  defect 
formation  in  an  inverse  transition  (W  to  S). 

In  order  to  study  the  spatial  variability  of  this  impulse  regime,  we  visualize  the  spatial  evolution  of  the  degree 
k{  of  the  zth  node,  using  polar  coordinates  (rf,  N  d  where  rf  =  k{  and  indicate  the  position  of  the  node 

which  is  fixed  on  the  outer  circumference  of  the  cylindrical  sample;  Fig.  3c.  We  refer  to  these  configurations  as 
“K-strings”  and  the  normal  vector  of  the  K- strings  at  each  node  indicates  the  local  direction  of  increasing  or 
decreasing  k{  with  time.  We  evaluate  the  variation  of  r,-  =  kt  at  each  position  (node)  while  we  consider  the  temporal 
evolution  of  each  single  event  (Fig.  3;  Figs  S5  and  S6).  In  Fig.  3,  we  show  that  the  onset  of  the  impulse  zone  coin¬ 
cides  with  the  folding  of  K-strings  where  the  normal  vectors  are  flipped  at  the  onset  of  the  non-linear  regime  and 
form  a  local  domain  that  we  refer  to  as  a  “kink”  (Fig.  3d).  Formation  of  domains  in  the  course  of  the  S-W  transi¬ 
tion  results  in  a  non-linear  behaviour  in  the  S  phase  of  the  ^-profile.  Since  the  trend  of  R(t)  mirrors  the  mean 
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Figure  3.  Nucleation  of  kinks  and  formation  of  domains,  (a)  The  mean  number  of  edges  (<k>)  versus  time 
for  event  #24  from  Lab.EQ4.  (b)  Transition  to  nucleation  zone  is  imprinted  in  diverging  the  inverse  of  mean 
betweeness  centrality  (B.C).  (c)  Schematic  representation  of  sensors  location  (red  filled  circles)  where  the  radius 
of  the  ring  is  proportional  with  node’s  degree  (d)  We  have  shown  accumulated  2D  spatio-temporal  patterns  of 
nodes’  degree  in  the  polar  system  for  each  time  interval  as  in  panel  (a).  Transition  from  the  linear  stage  (1)  to 
the  non-linear  regime  (2  +  3)  is  indicated  by  the  onset  of  local  defects  (black  arrows),  inducing  formation  of 
pair-kinks.  In  the  example  reported,  there  are  four  major  defect-zones.  The  arrows  are  normal  to  strings  and 
crumpled  strings  destroy  long-range  order  in  string  normal.  We  have  shown  more  examples  in  Figs  S5  and  S6. 


value  of  k  for  all  nodes  <k(t)>  (Fig.  S4),  we  might  also  infer  <k  (t)  >  oc  a  (t).  Given  this  observation,  it  is  clear 
that  the  non-linearity  of  stress  is  well-connected  to  nucleation  of  these  network  kinks. 

As  emphasized  by  Polyakov14  in  the  context  of  string  theory,  crumpled  stings  are  analogous  with  the 
Heisenberg  paramagnet  while  undulations  destroy  long-range  order  in  surface  normals14,15.  These  topological 
defects  are  local  defects  in  initially  ordered  structures  and  can  be  removed  by  global  collapse  of  K-strings  and 
local  bending  or  twisting  around  the  defects  cannot  remove  them  (i.e.,  they  are  topological  defects)16,17.  To  ana¬ 
lyze  deformation  of  K-strings,  we  map  them  onto  simplified  spin-like  chains  where  for  each  node  we  assign 
sf  =  sign(^pj  and  then  s2  =  ±1.  With  this  mapping,  defects  represented  with  negative  s2  indicate  flipped, 

inward-pointing  normal  vectors  (spins).  We  have  simplified  a  true  3D  configuration  of  acoustic  networks  by 
assigning  one  component  per  each  node  (up  or  down)  -see  Fig.  S7.  A  double  kink  separating  zones  with  up  and 
down  spins  functions  as  a  locator  for  the  change  from  one  ground  state  (S-phase)  to  another  degenerate  ground 
state  (W-phase)17-19. 

The  critical  point  < k>max  is  defined  when  \m\  =  <s>  approaches  its  minimum  value  (Figs  4  and  5b).  Here,  m 
is  the  order  parameter  of  the  K-strings  and  the  transition  from  S  to  W  (and  vice  versa)  occurs  continuously  A 
stable  disorder  phase  m  «  0  precedes  the  onset  of  the  S-phase,  and  the  system  of  nodes  is  forced  from  a  disordered 
state  (prior  to  S-phase)  to  an  ordered  state  (S).  This  happens  continuously  and  is  a  symmetry-breaking  transition. 
The  disordered  state  is  a  symmetric  state  and  in  the  ordered  phase  the  order  parameter  (m)  chooses  one  direction; 
for  the  S-phase  the  mean  direction  is  positive  ( j  outward),  whereas  for  the  W-phase,  the  mean  direction  is  nega¬ 
tive  ( l  inward).  Approaching  </c>max,  the  system  rapidly  transverses  a  temp  or  ary- unstable  symmetric  state.  The 
reason  lies  in  the  fact  that  the  symmetric  state  is  not  the  state  of  minimum  energy  and  that  in  the  process  of  evolv¬ 
ing  toward  the  ground  state,  the  symmetry  of  the  system  has  been  broken18. 

To  proceed,  we  calculate  a  correlation  function  for  the  K-strings  as  they  approach  < /c>max  (Fig.  4a, b).  We  can 
fit  a  correlation  function  such  as  G(x)  =  ^1  —  y )  exp  where  L  is  the  total  number  of  nodes,  x  is  distance, 
and  £  is  the  correlation  length19,20.  Correlation  length  £  is  tne  cut-offlength  of  where  for  distances  shorter  than 


SCIENTIFIC  REPORTS  |  6:21210  |  DOI:  10.1038/srep21210 


4 


www.nature.com/scientificreports/ 


a) 


b) 


s 

K 


Figure  4.  Continuous  Phase  Transition  from  S  to  W.  Results  from  a  cracking  noise  #255  from  Lab.EQ4  with 
mapping  on  the  Ising-chain  (a)  The  real-time  order  parameter  m  versus  time,  (b)  Real-time  cross-correlation 
function.  Very  close  to  the  transition  point  in  red,  we  can  see  nearly  overlapping  patterns  of  G(x)  indicating  a 
frozen  zone  of  correlation  length  (c)  cross- correlation  function  versus  the  normalized  distance.  For  fully 
ordered  state  a  triangular  function  (fully  coherent)  is  given  by  green  line.  Approaching  transition  point,  we  can 
fit  a  correlation  function  (red  line),  to  determine  correlation  length.  For  this  event,  the  frozen  correlation  length 
is  |  ^  0.086.  Blue  points  are  the  experimental  measures  of  the  correlation  function.  Before  normalization,  we 
had  F  =  300  nodes  and  |  «  21  nodes.  Approaching  the  critical  point,  the  correlation  length  becomes  frozen  as 
shown  in  (d).  For  this  event,  the  correlation  length  is  roughly  constant  for  —0.6  pus.  See  supplementary  Figures 
S13-S15  for  more  examples. 


f  the  correlation  function  G(x)  can  be  fitted.  For  a  fully  ordered  state  a  triangular  function  (where  f  — ►  «>)  is  given 
by  the  green  line  in  Fig.  4c. 

Interestingly  the  correlation  length  becomes  essentially  frozen  as  </c>max  is  approached,  as  shown  in  Fig.  4d 
for  crackling  event  #255.  More  examples  shown  in  Fig.  5  indicate  that  the  results  are  universal  for  the  recorded 
acoustic  events.  This  means  as  the  system  approaches  the  critical  point  at  a  finite  rate,  and  after  some  point  in  time 
the  correlation  length  cannot  maintain  its  equilibrium  value  and  a  transition  occurs  with  f  much  smaller  than 
the  system  size.  This  correlation  length  sets  the  mean  finite-size  of  the  final  domains.  The  formation  of  domains 
as  well  as  the  observation  of  frozen  correlation  length  at  the  critical  point  of  our  acoustic  networks  indicate  the 
existence  of  out- of- equilibrium  mechanism  that  is  well-described  by  the  Kibble-Zurek  mechanism  (KZM,  see 
Methods  section)7-10. 

The  core  idea  behind  the  KZM  is  that  near  the  phase  transition,  freezing  of  the  correlation  length  is  unavoid¬ 
able.  Based  on  this  theory,  the  resulting  density  of  defects  left  behind  by  continuous  transitions  is  dependent  on 
the  rate  at  which  the  critical  point  is  traversed,  and  the  rate  with  which  the  system  can  adjust,  defining  the  relaxa¬ 
tion  or  healing  time  of  the  system.  This  mechanism  is  reflected  in  the  density  of  defects  and  the  “ freeze-out ”  time 
which  scales  with  ramp -rate.  To  test  the  KZM  s  defect  density  prediction,  we  carefully  measured  the  number  of 
flipped  nodes  in  the  vicinity  of  the  critical  point  where  correlation  length  is  frozen  (Fig.  6). 

A  key  output  of  our  analysis  is  that  the  number  of  defects  (i.e.,  flipped  nodes  at  the  final  state)  is  larger  when 

the  local  ramp  rate  ^  is  faster  (Fig.  6a).  We  can  fit  a  power-law  scaling  as:  |  oc  ^  j  (Fig.  6a),  where  |  is  the 

frozen  correlation  length  and  ( — )  is  the  ramp  rate  (Figs  S4-6  and  S10).  To  determine  the  ramp  rate,  which  is 

V  dt  )t=tc 

analogous  to  the  local  loading  rate  prior  to  the  nucleation  of  kinks,  we  measure  the  slope  of  R(t)  (Fig.  3b).  The 
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Figure  5.  Real-time  correlation  length  of  Acoustic-Networks.  (a,b)  The  real-time  evolution  of  the  order 
parameter  in  the  S-phase  and  the  corresponding  pair- correlation  for  a  dry-cracking  event  from  Lab.Eq3.  The 
fully  coherent  system  has  a  triangular  shape.  rnorm  is  the  normalized  distance  between  nodes.  For  early  time 
steps  (t  <  3  |is),  thermal  activation  of  kinks  does  not  destroy  the  long-range  correlation  which  spans  the  whole 
system  (c,d)  When  the  system  approaches  the  critical  point,  the  correlation  length  £  is  effectively  “frozen”  as 
shown  by  the  roughly  horizontal  portions  of  the  curves  in  C.  The  size  of  the  network  system  was  300  nodes  and 
for  the  acoustic-networks  far-beyond  the  critical  point  a  long-range  correlation  length  is  assigned,  i.e.,  the 
coherence  spans  the  whole  system. 
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Figure  6.  Dependency  of  frozen  correlation  length  (kink  density)  on  ramping  rate,  (a)  Events  with  faster 
transition  to  their  critical  points  induce  higher  defect  density  (i.e.,  shorter  correlation  length).  Here  we  show 
typical  rupture  fronts  from  Lab.EQ3  and  the  size  of  the  network  is  300  nodes.  We  obtained  b  «  0.35  in 
f  oc  (Rc)  (dashed  line),  in  agreement  with  the  mean-field  model  prediction  (also  see  Figs  S14-S16)  (b)  effect 
of  local  stress-ramp  rates  on  the  fast- weakening  regime:  Events  with  faster  weakening  rate  scale  with  slower 
local  ramp  rate.  This  is  illustrated  here  via  a  log-log  plot  of  the  normalized  rate  of  weakening  Rw  observed  as  a 
function  of  ramp-rate  (events  from  Lab.EQ4). 
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exponent  of  ~0.35  ±  0.06  obtained  by  fitting  the  laboratory  data  is  in  agreement  (within  experimental  error)  with 
the  theoretical  value  v/(l  +  vz)  =  0.34  where,  for  the  2d  classical  Ising  system,  v=  1  and  z=  2. 119,21.  Here,  the 
parameters  u  and  z  are  spatial  and  dynamical  critical  exponents  (see  Methods  section). 

Furthermore,  we  can  evaluate  the  time-reversal  transition  from  the  W  to  S  phase  while  we  approach  to  the 
critical  point  from  right  (in  other  words,  if  we  heal  or  reverse  the  failure  process).  Approaching  from  the  left 
(S-phase)  or  right  (W-phase;  time  reversal  or  healing  scenario)  to  the  critical  point  Ac  results  in  slightly  different 
characteristics  of  the  defect  density  (Fig.  5b).  The  rate  of  the  S-ramp  (linear  strengthening  rate  or  ramp  rate)  for 
most  of  the  recorded  events  is  higher  than  the  rate  of  W-ramp  (i.e.,  linear  weakening  rate).  Next  we  estimate  the 
linear  weakening  rate.  Using  the  KZM  scaling  law  for  defects,  we  obtain  (Supplementary  Information): 

1  +UZ 

—2v 

,  in  which  r0  is  determined  by  the  microscopic  details  of  the  system;  rs  and  rB,  are  the 


1 


(T 


-2/7 
+  UZ 


time  characteristics  of  ramps  in  the  S  and  W  phases,  respectively  Therefore,  the  linear  weakening  rate  is  given  by: 
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.  A  faster  local  S-ramp  inversely  scales  with  weakening  rate.  This  is  an  important 


result  since  it  has  been  shown  that  weakening  rate  is  correlated  with  the  global  rupture  velocity  of  cracks3,12,22.  To 
verify  this  prediction,  we  measured  the  rate  of  weakening  from  R  profiles  R  =  Rmax  ~  Rmin .  As  shown  in  Fig.  6b, 

W  Tw 

our  measurements  confirm  the  aforementioned  prediction. 


Discussion 

We  have  shown,  for  the  first  time  via  laboratory  data,  the  evolution  of  propagating  microcracks  over  the  duration 
of  only  a  few  microseconds  derived  from  multiple  AE  event  data.  We  mapped  acoustic  excitations  from  crackling 
events  to  complex  networks  and  then  further  to  spin-like  chains.  Applying  these  tools  in  a  novel  new  way  to  AE 
data,  we  elucidated  the  transition  to  from  a  precursory  strengthening  phase  to  a  phase  of  rapid  weakening  and 
unstable  crack  growth.  We  demonstrated  a  real-time  probe  of  evolving  ground  state  of  the  system  via  acoustic 
(phononic)  excitations  and  then  we  accessed  the  strongly  non-equilibrium  dynamics  directly,  rather  through  its 
aftermath.  We  illustrated  that  as  we  force  our  system  to  cross  the  transition  at  a  finite  velocity,  different  regions 
of  the  system  will  choose  different  minima  of  the  free  energy  This  leads  to  the  appearance  of  topological  defects. 
Using  an  analogy  with  the  Ising  ferrormagnet  model,  we  discovered  that  topological  defects  correspond  to  the 
appearance  of  regions  where  the  nodes  form  domains  pointing  either  up  or  down. 

As  we  described,  flipped  nodes  induce  a  degree  of  non-linearity  (inelastic  behaviour)  in  dynamic  stress  pro¬ 
files,  and  dislocations  are  good  physical  candidates  for  such  defects  in  crystalline  solids.  Far  from  the  critical  point 
where  the  system  is  still  coherent,  thermal  activation  of  dislocations  occurs  adiabatically;  however  rapid-traversing 
to  the  failure  point  violates  this  assumption.  As  another  perspective  to  support  nonlinearity  of  acoustic  emission 
waves,  we  can  assume  that  a  K- string  is  a  classical  vibrating  string  (i.e.,  harmonic  or  anharmonic  oscillators) 
where  “folding”  adds  normal  modes  to  the  oscillator23,24.  Assuming  a  string  with  the  fixed  boundaries  at  the  end 
of  elastic  string  has  n  normal  modes,  the  boundary  value  problem  of  solving  the  wave  equation  results  in: 
un  (x,  t)  =  exp  ( —  iuj  Kt)  sin  — x,  where  ujk  is  the  frequency  of  an  excitation  of  wave- vector  K,  L  is  the  fixed  length 
of  the  string  (number  of  nocles),  and  un(x,  t )  represents  the  deflection  of  the  string  which  satisfies  the  boundary 
condition.  The  fundamental  mode  n=  1  for  our  network- strings  is  a  fully-ordered  state.  Each  flipped  node  is  an 
analogy  with  phononic  (or  generally  bosonic)  excitation.  We  might  call  them  “ nodons ”  to  specify  excitations  on 
network-like  structures.  For  a  given  wave-vector  ( K ),  the  string  could  have  nK  nodons  (Fig.  S12).  With  this 
description,  we  support  the  previous  speculation  on  the  non-linear  nature  of  acoustic  waves  emitted  from  differ¬ 
ent  sources25,26. 

While  our  primary  focus  in  this  study  was  to  characterize  events  with  the  definite  continuous  S-W  transition, 
we  have  also  recognized  events  with  an  abrupt  change  in  the  order  parameter  that  is  characteristic  of  a  first  order 
transition  (Fig.  S18).  Further  study  is  needed  to  explore  more  features  of  first-order  “laboratory  earthquakes”. 
In  addition,  analysis  of  3D  acoustic  networks  by  mapping  them  to  spin  models  will  present  a  unique  picture  of 
the  evolution  of  microcracks  under  true  3D  stress-fields.  It  would  be  interesting  to  study  the  transition  from 
fast-weakening  phase  to  the  next  phase  and  monitor  the  evolution  of  nodes’  states  under  a  true  3D  acoustic 
networks  (such  as12)  where  this  transition  defines  the  crack-like  or  pulse-like  nature  of  rupturing27.  Another 
interesting  study  could  focus  on  the  dynamics  of  annihilation  of  flipped  nodes  and  its  effect  on  self-healing  pulses. 
Finally,  our  results  could  be  extended  in  the  study  of  dynamics  of  frictional  interfaces  where  dissipation  of  energy 
is  coupled  with  the  variation  of  contact  areas13  as  well  as  the  study  of  stick-slip  motion  with  the  formation  of 
kinks  and  antikinks26. 


Methods 

Laboratory  Procedures.  We  use  four  sets  of  recorded  acoustic  emissions  (labeled  as  Lab.EQl,  2,  3  and  4) 
from  Westerly  granite  and  Basalt  rock  samples  (most  of  the  analyzed  events  are  precursor  rupture  fronts).  The 
Lab.EQl  and  2  are  the  recorded  multi- stationary  acoustic  waveforms  from  evolution  of  frictional  rock- interfaces 
of  Westerly  Granite  samples.  The  interfaces  were  in  dry  conditions  with  smooth  (saw-cut)  and  naturally  rough 
surfaces,  respectively28.  The  evaluated  events  are  from  different  stages  and  position  and  are  not  limited  to  par¬ 
ticular  stage  of  the  tests.  The  Lab.EQ3  is  the  fast-loading  experiment  on  a  cylindrical  sample  of  Westerly  Granite 
(~  10-5  s_1)  at  50  MPa  confining  pressure  (approximately  2  km),  which  is  about  an  order  faster  than  Lab.EQl  and 
229  Lab.EQ4  are  events  from  Basalt  samples;  described  in30.  The  global  loading  rate  was  10_6s_1.  In  all  of  the  above 
experiments,  we  reordered  amplified  events  using  16  to  18  sensor  networks  in  both  short  (discrete  events)  and 
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long  timescale  recorders  (AE  records).  The  resolution  of  each  recorded  interval  during  the  life-time  of  a  waveform 
was  -20-100  ns. 

Networks  of  Acoustic  Emission  Waveforms  .  The  concept  behind  studying  each  single  acoustic  excita¬ 
tion  event  -  in  this  study  -  is  to  characterize  sub-events  involved  in  the  course  of  just  a  single  AE  event.  A  big 
avalanche  is  composed  of  many  smaller  components,  which  trigger  one  another1.  To  study  each  acoustic  event, 
we  use  functional  network  theory  to  analyze  multiple  recorded  waveforms. 

To  evaluate  reordered  multiple  acoustic  emissions  (multiple  time  series  for  a  single  event),  we  use  a  previous 
algorithm  on  waveforms  from  our  reordered  acoustic  emissions4,11,12.  The  main  steps  of  the  algorithm  are  as 
follows11: 

(1)  The  waveforms  recorded  at  each  acoustic  sensor  are  normalized  by  the  maximum  value  of  the  amplitude  in 
that  station. 

(2)  Each  time  series  is  divided  according  to  maximum  segmentation,  in  a  way  that  each  segment  includes  only 
one  data  point.  The  amplitude  of  the  jth  segment  from  the  zth  time  series  (1  <  i  <  N)  is  denoted  by  u^{t)  (in 
mV).  N  is  the  number  of  nodes  or  acoustic  sensors. 

(3)  ul,i(t )  is  compared  with  uk^(t)  to  create  links  between  the  nodes  using  the  following  methodology:  If 

d  (zA;  (t) ,  uk^  (t) )  <  £  (where  £  is  the  threshold  level  discussed  in  the  following  point)  we  set  aik  (J)  =  1,  oth¬ 
erwise  aik(j)  =  0  where  aik(j )  is  the  component  of  the  connectivity  matrix  and  d  =  \\u1^  ( t ) ,  ( t )  ||  is  the 

employed  similarity  metric.  The  employed  norm  in  our  algorithm  is  the  absolute-norm.  With  this  metric,  we 
simply  compare  the  amplitude  of  sensors  at  each  time- step. 

(4)  Threshold  level  (£):  To  select  a  threshold  level,  we  use  a  method4,11  that  uses  an  adaptive  threshold  criterion 
and  is  stable  for  a  range  of  deviations  from  £.  The  result  of  this  algorithm  is  an  adjacency  matrix  with  compo¬ 
nents  given  by  a  (xz(t),  xk{t))  =  0(£  —  |  u1^  (t)  —  uk^  (t)  |)  where  O  (...)  is  the  Heaviside  function. 

In  general,  the  modularity  of  a  network  measures  the  degree  of  division  of  that  network  into  modules  (clus¬ 
ters):  if  a  network  has  high  modularity,  the  strength  of  connections  in  individual  modules  is  strong,  whereas  the 
strength  of  connections  between  modules  is  not.  The  networks  modularity  characteristic  is  addressed  as  the 
quantity  of  densely  connected  nodes  relative  to  a  null  (random)  model31.  The  modularity  is  quantified  using  the 
Q-profile ,  and  is  the  result  of  some  optimization  of  the  cluster  structure  of  a  given  network.  The  modularity  Q  (i.e., 
the  objective  function)  is  defined  as23: 


nm 

Q  =  E 

s= 1 


h_ 

Af 

L 

PL) 

(A.1) 


in  which  NM  is  the  number  of  modules,  L  =  kP  h  is  the  number  of  links  in  module  s  and  ds  =  JT  k-  (the 

sum  of  node  degrees  in  module  s). 

Then,  in  each  time  step  during  the  evolution  of  the  waveforms  (here  over  observation  windows  of  -400  |is),  we 
obtain  a  Q  value.  The  temporal  evolution  of  Q  in  the  monitored  time  interval  forms  the  Q-profile. 

This  network  algorithm  can  be  explained  in  the  context  of  space- correlation  methods.  We  can  define  a  similar 
measure  to  a  time- windowed  correlation  method32  where  the  inner  product  is  replaced  with  a  Heaviside  function. 
Let  us  consider  a  sequence  of  nodes  over  a  certain  time  step.  The  space-windowed  correlation  is  given  by32: 


«(*,) 


py  u{x')u{x'  +  xs)dx' 


rx+l 

Jx-l 


u  ( xr )  dxr 


(A.2) 


where  the  space- window  of  length  L=  2/  is  centered  at  length  /,  xs  is  the  space  shift  used  in  the  cross  correlation, 
and  the  amplitude  of  each  node  is  u(x).  Here  the  employed  norm  is  the  inner-product  and  can  be  replaced  with 
another  norm: 


n  x-\-l 

R{xs)  oc  J  \\u(x'),  u(x'  +  xs)\\dx'.  (Ai3) 

Summing  over  all  space- shifts  and  replacing  the  norm  with  our  similarity  metric  we  get: 

Toe EE0(C-  IK*'),  «(*'  +  xs)\\)  , 

xs  x'  (A.  4) 

T  is  proportional  to  density  of  links  determined  by  the  similarity  metric  used  to  construct  links  between  nodes 
represented  by  pairs  of  amplitudes  u  (xr) ,  u  (xf  +  xs). 

Kibble-Zurek  Mechanism.  The  idea  behind  the  KZM  is  to  compare  the  relaxation  time  (or  healing  time 
of  the  system  in  equilibrium)  with  the  timescale  of  change  in  the  control  parameter  (e).  Assuming  a  linear 
change  of  control  parameter  in  the  vicinity  of  the  critical  point  z(t)  =  t/rs ,  where  rs  is  the  ramp  time  in  S-phase. 

The  relaxation  or  healing  time  we  consider  is  an  equilibrium  (quasi-static)  condition:  r  (£)  =  and  vz=  p. 
This  determines  the  reaction  time  of  the  order  parameter.  Here,  v  and  z  are  spatial  and  dynamical  critical  expo¬ 
nents,  and  t0  is  a  characteristic  timescale8,9.  The  system  can  adiabatically  follow  the  change  imposed  by  the 
local  stress  ramp  if  the  relaxation  time  characterized  by  r(e)  is  outside  the  interval  set  by  the  “freeze-out”  time 
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t  =  (r0rsn  1+  vz  centered  around  the  transition  point  (see  Fig.  S3).  The  correlation  length  will  effectively 
remain  fixed  (i.e.,  “freeze”)  at  time  t  before  reaching  the  critical  point.  The  correlation  length  is  given  by8,9,20: 
|  =  £0(ts/t0)i+vz,  where  £0  is  a  characteristic  correlation  length. 

Topological  defects  are  formed  with  the  density  of  one  defect  fragment  per  domain.  An  estimate  for  the  result¬ 
ing  density  of  topological  defects  is  given  by8,9:  p  oc  (t0/ts)T+vz*  In  the  frozen  phase,  one  can  define  an  effec¬ 


tive  control  parameter  £  =  R  =  —  =  l  — 


(?) 


....  «. 

1/1+/I  ^  ^  1/1  + 

21,33.  By  plotting  events  in  R  —  Rc  space  in  which  R  oc  (Rc) 


we  obtain  vz=  ft.  Then,  we  can  estimate  v  from  f  oc  (. Rcy+vz  where  we  measure  the  frozen  correlation  length  for 
the  given  event  with  the  rate  of  transition  Rc.  This  procedure  leads  to  estimates  of  1  and  z=  2.1  which  agrees 
with  the  mean-field  approximation  of  the  scaling  coefficients  of  the  2D  Ising-model  (Z2  symmetry  breaking19,21). 
To  verify  the  KZM  prediction  for  the  scaling  exponent,  we  analyzed  many  events  for  which  the  order  parameter 
changes  continuously  (Fig.  S15-17). 
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